Auckland Islands SRW satellite tracks

1 Filtering satellite tracks

1.1 Load and format

# read all .csv of Locations and convert to spatial object
ff <- list.files("../Data/NZ_SRW_Argos_data/Download_13April22/All_NZ_locations/", full.names = T)

# convert to sf object
raw_tracks <- ldply(ff, read.csv) %>% 
  st_as_sf(coords = c("Longitude", "Latitude"), crs = 4326)

# reintegrate lat lon in columns because we will need them later
raw_tracks[c("Longitude", "Latitude")] <- st_coordinates(raw_tracks)

# project spatial points to an antarctic coordinate system
raw_tracks <- st_transform(raw_tracks, proj_antarctic) 

1.2 Clean data

Remove satellite locations with quality Z and apply speed filter

# how many raw locations do we have
raw_tracks %>% 
  group_by(DeployID) %>% 
  dplyr::summarize(nb_locations = n())

## Simple feature collection with 18 features and 2 fields
## Geometry type: MULTIPOINT
## Dimension:     XY
## Bounding box:  xmin: -1990937 ymin: -5204887 xmax: 6114844 ymax: 966054.4
## CRS:           +proj=stere +lat_0=-90 +lat_ts=-71 +lon_0=0 +k=1 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0
## # A tibble: 18 × 3
##    DeployID    nb_locations                                             geometry
##    <chr>              <int>                                     <MULTIPOINT [m]>
##  1 212499              1751 ((1042776 -4331716), (1042878 -4331340), (1042958 -…
##  2 212500               101 ((522744.7 -4335400), (527629.6 -4302724), (527647.…
##  3 215258              2131 ((1042341 -4321852), (1047942 -4337869), (1048065 -…
##  4 215259              1046 ((-191707.8 -4113366), (930421.1 -4271630), (930693…
##  5 215260                 2            ((390322.3 -4915529), (2000106 -5204887))
##  6 215261              1170 ((1019345 -4226342), (1023114 -4316903), (1023437 -…
##  7 215262               522 ((687454.5 -4148217), (702342.2 -4128871), (705314.…
##  8 215263              2267 ((913058.8 -4310695), (919566.4 -4399070), (949488 …
##  9 46633               2670 ((934563.9 -4652730), (939939.8 -4656590), (940692.…
## 10 46635               1204 ((1034519 -4391005), (1035491 -4393757), (1036133 -…
## 11 46950               3983 ((-1422054 -4268807), (-1418681 -4263496), (-141107…
## 12 46955               1797 ((353001.1 -3093163), (354964.3 -3149689), (355070.…
## 13 Bill Wiremu         8359 ((-1990937 -2646971), (-1111275 -5027211), (814011.…
## 14 Rima                1850 ((1008044 -4357885), (1008238 -4366671), (1008510 -…
## 15 Rua                  243 ((1012690 -4362642), (1012977 -4365951), (1017595 -…
## 16 Tahi                2969 ((-1892230 -5009731), (837729.4 -4381394), (890273.…
## 17 Toru                3560 ((1025301 -4338445), (1026972 -4339054), (1027273 -…
## 18 Whaa                1660 ((753046.5 -4163074), (761742.7 -3922445), (761757.…

# tag 215260 only has two locations so I am removing it
raw_tracks <- raw_tracks %>% filter(DeployID != "215260")

# filter z positions
raw_tracks <- raw_tracks %>% filter(Quality != 'Z')

# transform Date into POSIX format
raw_tracks <- raw_tracks %>% 
  mutate(Date = as.POSIXct(Date, format = "%H:%M:%S %d-%b-%Y"),
        month = as.numeric(strftime(Date, format = "%m")))

# prefilter positions with the speed filter (see Reisinger et al. 2021 who did that before actually fitting the ssm)
# vmax in meters / s
raw_tracks <- ddply(raw_tracks, ~DeployID, function(d){
    d$argosfilter <- sdafilter(lat = d$Latitude, 
                               lon = d$Longitude, 
                               lc = d$Quality, 
                               dtime = d$Date, vmax = 5)
    return(d)
})
raw_tracks <- raw_tracks %>% 
  filter(argosfilter != "removed") %>% 
  dplyr::select(-argosfilter)

# who is left? how many positions?
raw_tracks %>% 
  group_by(DeployID) %>% 
  dplyr::summarize(nb_locations = n())

## # A tibble: 17 × 2
##    DeployID    nb_locations
##    <chr>              <int>
##  1 212499              1550
##  2 212500                95
##  3 215258              1829
##  4 215259               903
##  5 215261              1008
##  6 215262               430
##  7 215263              1980
##  8 46633               2413
##  9 46635               1053
## 10 46950               3486
## 11 46955               1679
## 12 Bill Wiremu         7330
## 13 Rima                1443
## 14 Rua                  207
## 15 Tahi                2507
## 16 Toru                2989
## 17 Whaa                1472

g <- gbase + 
  geom_sf(data = raw_tracks, aes(col = DeployID, geometry = geometry), size = 0.8)
ggsave(g, file = "./Outputs/filtered_tracks_NZ.png", dpi = 300, width = 8, height = 8)  
g

Figure 1.1: Filtered satellite tracks of right whales tagged in the Auckland Islands in 2020 and 2021.

1.3 Segmentate when transmissions interrupted

Assess mean time step between locations

timelaps_df <- ddply(raw_tracks, ~DeployID, function(d){
  d$timelaps_hours <- NA
  for (i in 2:nrow(d)){
    d$timelaps_hours[i] = difftime(d$Date[i], d$Date[i-1], units = "hours")}
  return(d)
})

# mean time step between locations (in hours)
a <- aggregate(timelaps_hours~DeployID, timelaps_df, mean)
a

##       DeployID timelaps_hours
## 1       212499      1.6367726
## 2       212500      4.8701537
## 3       215258      2.9402048
## 4       215259      2.8828732
## 5       215261      1.7521861
## 6       215262     15.0779429
## 7       215263      2.0706410
## 8        46633      1.5360574
## 9        46635      1.9188844
## 10       46950      1.8284794
## 11       46955      0.9882587
## 12 Bill Wiremu      1.2479754
## 13        Rima      1.2587930
## 14         Rua      4.6188660
## 15        Tahi      3.0898524
## 16        Toru      1.3542783
## 17        Whaa      0.9798064

mean(a$timelaps_hours)

## [1] 2.944237

ggplot(timelaps_df, aes(timelaps_hours)) + 
  geom_histogram(binwidth = 5, col ="white", na.rm = T) + 
  theme_bw() + 
  facet_wrap(~DeployID, scales = "free") +
  xlab("Time laps between successive locations (in hours)")

Figure 1.2: Interval between positions recorded by the tags (in hours). These histograms are used to assess potential gaps in the transmission and to set the time step in the interpolated tracks.

source("Fun_TrackInterruption_SRW.R")
# need a time column to run
raw_tracks$time <- raw_tracks$Date
# run the custom function to split tarcks interrupted for more than 144 hours = 6 days
raw_tracks <- Fun_TrackInterruption(dataframe = raw_tracks, hours_gap = 144, id = "DeployID")

# set first timelaps of each segment to NA
raw_tracks <- ddply(raw_tracks, ~segmentid, function(d){
  d$timelaps[1] <- NA
  return(d)
})

# print number of locations per segments
raw_tracks %>% 
  group_by(segmentid) %>% 
  dplyr::summarize(nb_locations = n()) %>% 
  print(n = Inf)

## # A tibble: 27 × 2
##    segmentid   nb_locations
##    <chr>              <int>
##  1 212499              1550
##  2 212500_a              84
##  3 212500_b              11
##  4 215258_a            1727
##  5 215258_b             102
##  6 215259_a             889
##  7 215259_b              14
##  8 215261_a             971
##  9 215261_b              37
## 10 215262_a             376
## 11 215262_b              30
## 12 215262_c               7
## 13 215262_d              17
## 14 215263              1980
## 15 46633               2413
## 16 46635               1053
## 17 46950               3486
## 18 46955               1679
## 19 Bill Wiremu         7330
## 20 Rima                1443
## 21 Rua_a                196
## 22 Rua_b                 11
## 23 Tahi_a              2195
## 24 Tahi_b                 2
## 25 Tahi_c               310
## 26 Toru                2989
## 27 Whaa                1472

# remove segments with less than 15 locations
toremove <- raw_tracks %>% 
  group_by(segmentid) %>% 
  dplyr::summarize(nb_locations = n()) %>% 
  filter(nb_locations <= 15)

raw_tracks <- raw_tracks %>% filter(!(segmentid %in% toremove$segmentid))

# print number of locations per segments
raw_tracks %>% 
  group_by(segmentid) %>% 
  dplyr::summarize(nb_locations = n()) %>% 
  print(n = Inf)

## # A tibble: 22 × 2
##    segmentid   nb_locations
##    <chr>              <int>
##  1 212499              1550
##  2 212500_a              84
##  3 215258_a            1727
##  4 215258_b             102
##  5 215259_a             889
##  6 215261_a             971
##  7 215261_b              37
##  8 215262_a             376
##  9 215262_b              30
## 10 215262_d              17
## 11 215263              1980
## 12 46633               2413
## 13 46635               1053
## 14 46950               3486
## 15 46955               1679
## 16 Bill Wiremu         7330
## 17 Rima                1443
## 18 Rua_a                196
## 19 Tahi_a              2195
## 20 Tahi_c               310
## 21 Toru                2989
## 22 Whaa                1472

# now what's the average time between locations?
a <- aggregate(timelaps~segmentid, raw_tracks, mean)
a

##      segmentid  timelaps
## 1       212499 1.6367726
## 2     212500_a 2.2604284
## 3     215258_a 2.2659622
## 4     215258_b 2.8483443
## 5     215259_a 2.6024343
## 6     215261_a 1.4457380
## 7     215261_b 5.1853704
## 8     215262_a 2.3169081
## 9     215262_b 5.2023563
## 10    215262_d 9.8943056
## 11      215263 2.0706410
## 12       46633 1.5360574
## 13       46635 1.9188844
## 14       46950 1.8284794
## 15       46955 0.9882587
## 16 Bill Wiremu 1.2479754
## 17        Rima 1.2587930
## 18       Rua_a 1.8735826
## 19      Tahi_a 2.4392563
## 20      Tahi_c 1.6405735
## 21        Toru 1.3542783
## 22        Whaa 0.9798064

mean(a$timelaps)

## [1] 2.490691

# what's the maximum time between locations?
a <- aggregate(timelaps~segmentid, raw_tracks, max)
a

##      segmentid timelaps
## 1       212499 22.39028
## 2     212500_a 21.43056
## 3     215258_a 70.58194
## 4     215258_b 36.63500
## 5     215259_a 49.24583
## 6     215261_a 22.82778
## 7     215261_b 37.23389
## 8     215262_a 34.74778
## 9     215262_b 21.46528
## 10    215262_d 51.34000
## 11      215263 22.48583
## 12       46633 21.35639
## 13       46635 42.47000
## 14       46950 18.19167
## 15       46955 10.59694
## 16 Bill Wiremu 23.38111
## 17        Rima 73.75361
## 18       Rua_a 20.35194
## 19      Tahi_a 22.30056
## 20      Tahi_c 15.07583
## 21        Toru 61.70472
## 22        Whaa 10.33222

save(raw_tracks, file = "./Outputs/raw_tracks.RData")

g <- gbase + 
  geom_sf(data = raw_tracks, aes(col = segmentid, geometry = geometry), size = 0.8)+
  theme(legend.position = "right")

ggsave(g, file = "./Outputs/filtered_tracks_NZ_segments.png", dpi = 300, width = 8, height = 8)
g

Figure 1.3: Filtered satellite tracks of right whales tagged in the Auckland Islands in 2020 and 2021 and South Georgia in 2020. Any segment of less than 15 positions separated by interruptions of more than 6 days was removed

2 State Space Model

2.1 Prepare data

Prepare data in a specific format for foieGras package

raw_tracks$time <- raw_tracks$Date

data_ssm <- raw_tracks %>% 
  dplyr::select(segmentid, time, Quality, Longitude, Latitude, Error.Semi.major.axis, Error.Semi.minor.axis, Error.Ellipse.orientation) %>% 
  dplyr::rename(id = segmentid, date = time, lc = Quality, lon = Longitude, lat = Latitude, smaj = Error.Semi.major.axis, smin = Error.Semi.minor.axis, eor = Error.Ellipse.orientation)
data_ssm$geometry <- NULL

str(data_ssm)

## 'data.frame':    32329 obs. of  8 variables:
##  $ id  : chr  "212499" "212499" "212499" "212499" ...
##  $ date: POSIXct, format: "2021-07-21 09:26:18" "2021-07-21 09:29:49" ...
##  $ lc  : chr  "A" "B" "A" "B" ...
##  $ lon : num  166 166 166 166 166 ...
##  $ lat : num  -50.5 -50.5 -50.5 -50.5 -50.5 ...
##  $ smaj: int  NA 1110 1751 6119 1359 1769 2128 4656 4081 5304 ...
##  $ smin: int  NA 222 43 117 802 1332 892 1109 1121 758 ...
##  $ eor : int  NA 61 106 172 166 93 117 93 98 107 ...

save(data_ssm, file = "./Outputs/data_ssm.RData")

2.2 Fit model with two different time steps

Fit spate-space model (based on codes from Reisinger et al. 2021)

nz_fit_ssm <- tibble(time_step = c(3, 6)) %>% # create a tibble with the two time step conditions in one column
  group_by(time_step) %>% # 
  nest() %>% 
  mutate(data = list(data_ssm, data_ssm), # store the same dataframe in the two rows
         fit = map2(time_step, data, function(x, y){
                     fit_ssm(y,
                        spdf = F, # we already run the filter so no need to do again
                        pf = F,
                        model = "rw",
                        min.dt = 48,
                        time.step = x,
                        control = ssm_control(verbose = 0),
                        map = list(psi = factor(NA)))}))

2.3 ARS

Assess movement persistence along tracks

nz_fit_ssm <- nz_fit_ssm %>% 
  mutate(fmp = map(fit, function(a){
    fit_mpm(a, what = "predicted", model = "mpm", control = mpm_control(verbose = 0))
  }))

save(nz_fit_ssm, file = "./Outputs/nz_fit_ssm.RData")

load("./Outputs/nz_fit_ssm.RData")

pdf(file = "./Outputs/nz_fit_ssm_interpolation.pdf", width = 12, height = 12)
nz_fit_ssm %>% 
  pull(fmp) %>% 
  llply(., function(x) {plot(x, pages = 1, ncol= 5, pal = "Zissou1", rev = T)})

## [[1]]

## 
## [[2]]

dev.off()

## png 
##   2

nz_fit_ssm %>% 
  pull(fmp) %>% 
  llply(., function(x) {plot(x, pages = 1, ncol= 5, pal = "Zissou1", rev = T)})

## [[1]]

Figure 2.1: Time series of movement persistence modeled along right whale tracks from Auckland Islands. The move persistence parameter gt ranges continuously from 0 (little persistence, indicative of area-restricted movements) to 1 (high persistence, indicative of directed movements).

## 
## [[2]]

Figure 2.2: Time series of movement persistence modeled along right whale tracks from Auckland Islands. The move persistence parameter gt ranges continuously from 0 (little persistence, indicative of area-restricted movements) to 1 (high persistence, indicative of directed movements).

3 Behavior State

3.1 Classify move persistence

Convert predictions to a dataframe and extract ARS

nz_fit_ssm_df <- nz_fit_ssm %>% 
  unnest(fmp) %>% 
  dplyr::select(-c(fit, data)) %>% 
  mutate(pred = map(mpm, function(x){
    cbind(x$data[c("date", "x", "y")], x$fitted["g"])
  })) %>% 
  dplyr::select(c(time_step, id, pred)) %>% 
  unnest(pred) %>% 
  st_as_sf(coords = c("x", "y"), crs = 4326) %>% 
  st_transform(proj_antarctic) 

# distribution of g values
summary(nz_fit_ssm_df$g)

##     Min.  1st Qu.   Median     Mean  3rd Qu.     Max. 
## 0.003713 0.418424 0.640095 0.587120 0.771512 0.976121

# select the 50% quantile as a threshold for ARS. it matches pretty well the cut-off used by Zerbini et al. with the beta parameter of the SSM
nz_fit_ssm_df <- nz_fit_ssm_df %>% 
    mutate(mode = case_when(g < quantile(g, 0.5) ~ "ARS",
                            g >= quantile(g, 0.5) ~ "transit"))

4 Compare time steps

fun_plot_ars <- function(ids) {
  for(i in ids){
    d <- nz_fit_ssm_df %>% filter(id == i)
    g1 <- ggplot(d, aes(col = g)) +
      geom_sf() +
      facet_wrap(~time_step) +
      scale_color_viridis_c(name = "Move persistence")
    g2 <- ggplot(d, aes(col = mode)) +
      geom_sf() +
      geom_sf(data = d[d$mode == "ARS", ], aes(col = mode)) + # plot again on top
      facet_wrap(~time_step)
    g <- (g1 / g2) &
      plot_annotation(title = i)
    print(g)
  }
}

pdf(file = "./Outputs/nz_fit_ssm_interpolation_tracks.pdf", width = 10, height = 8)
fun_plot_ars(unique(nz_fit_ssm_df$id))
dev.off()

## png 
##   2

fun_plot_ars(unique(nz_fit_ssm_df$id))

Figure 4.1: Movement persistence modeled along right whale tracks from Auckland Islands. The move persistence parameter gt ranges continuously from 0 (little persistence, indicative of area-restricted movements) to 1 (high persistence, indicative of directed movements).

Figure 4.2: Movement persistence modeled along right whale tracks from Auckland Islands. The move persistence parameter gt ranges continuously from 0 (little persistence, indicative of area-restricted movements) to 1 (high persistence, indicative of directed movements).

Figure 4.3: Movement persistence modeled along right whale tracks from Auckland Islands. The move persistence parameter gt ranges continuously from 0 (little persistence, indicative of area-restricted movements) to 1 (high persistence, indicative of directed movements).

Figure 4.4: Movement persistence modeled along right whale tracks from Auckland Islands. The move persistence parameter gt ranges continuously from 0 (little persistence, indicative of area-restricted movements) to 1 (high persistence, indicative of directed movements).

Figure 4.5: Movement persistence modeled along right whale tracks from Auckland Islands. The move persistence parameter gt ranges continuously from 0 (little persistence, indicative of area-restricted movements) to 1 (high persistence, indicative of directed movements).

Figure 4.6: Movement persistence modeled along right whale tracks from Auckland Islands. The move persistence parameter gt ranges continuously from 0 (little persistence, indicative of area-restricted movements) to 1 (high persistence, indicative of directed movements).

Figure 4.7: Movement persistence modeled along right whale tracks from Auckland Islands. The move persistence parameter gt ranges continuously from 0 (little persistence, indicative of area-restricted movements) to 1 (high persistence, indicative of directed movements).

Figure 4.8: Movement persistence modeled along right whale tracks from Auckland Islands. The move persistence parameter gt ranges continuously from 0 (little persistence, indicative of area-restricted movements) to 1 (high persistence, indicative of directed movements).

Figure 4.9: Movement persistence modeled along right whale tracks from Auckland Islands. The move persistence parameter gt ranges continuously from 0 (little persistence, indicative of area-restricted movements) to 1 (high persistence, indicative of directed movements).

Figure 4.10: Movement persistence modeled along right whale tracks from Auckland Islands. The move persistence parameter gt ranges continuously from 0 (little persistence, indicative of area-restricted movements) to 1 (high persistence, indicative of directed movements).

Figure 4.11: Movement persistence modeled along right whale tracks from Auckland Islands. The move persistence parameter gt ranges continuously from 0 (little persistence, indicative of area-restricted movements) to 1 (high persistence, indicative of directed movements).

Figure 4.12: Movement persistence modeled along right whale tracks from Auckland Islands. The move persistence parameter gt ranges continuously from 0 (little persistence, indicative of area-restricted movements) to 1 (high persistence, indicative of directed movements).

Figure 4.13: Movement persistence modeled along right whale tracks from Auckland Islands. The move persistence parameter gt ranges continuously from 0 (little persistence, indicative of area-restricted movements) to 1 (high persistence, indicative of directed movements).

Figure 4.14: Movement persistence modeled along right whale tracks from Auckland Islands. The move persistence parameter gt ranges continuously from 0 (little persistence, indicative of area-restricted movements) to 1 (high persistence, indicative of directed movements).

Figure 4.15: Movement persistence modeled along right whale tracks from Auckland Islands. The move persistence parameter gt ranges continuously from 0 (little persistence, indicative of area-restricted movements) to 1 (high persistence, indicative of directed movements).

Figure 4.16: Movement persistence modeled along right whale tracks from Auckland Islands. The move persistence parameter gt ranges continuously from 0 (little persistence, indicative of area-restricted movements) to 1 (high persistence, indicative of directed movements).

Figure 4.17: Movement persistence modeled along right whale tracks from Auckland Islands. The move persistence parameter gt ranges continuously from 0 (little persistence, indicative of area-restricted movements) to 1 (high persistence, indicative of directed movements).

Figure 4.18: Movement persistence modeled along right whale tracks from Auckland Islands. The move persistence parameter gt ranges continuously from 0 (little persistence, indicative of area-restricted movements) to 1 (high persistence, indicative of directed movements).

Figure 4.19: Movement persistence modeled along right whale tracks from Auckland Islands. The move persistence parameter gt ranges continuously from 0 (little persistence, indicative of area-restricted movements) to 1 (high persistence, indicative of directed movements).

Figure 4.20: Movement persistence modeled along right whale tracks from Auckland Islands. The move persistence parameter gt ranges continuously from 0 (little persistence, indicative of area-restricted movements) to 1 (high persistence, indicative of directed movements).

Figure 4.21: Movement persistence modeled along right whale tracks from Auckland Islands. The move persistence parameter gt ranges continuously from 0 (little persistence, indicative of area-restricted movements) to 1 (high persistence, indicative of directed movements).

Figure 4.22: Movement persistence modeled along right whale tracks from Auckland Islands. The move persistence parameter gt ranges continuously from 0 (little persistence, indicative of area-restricted movements) to 1 (high persistence, indicative of directed movements).

g <- gbase + 
  geom_sf(data = nz_fit_ssm_df, aes(col = g, geometry = geometry), size = 0.8) +
  theme(legend.position = "right") +
  scale_color_viridis_c(name = "Move persistence") +
  facet_wrap(~time_step)

ggsave(g, file = "./Outputs/nz_fit_ssm_tracks_all.png", dpi = 300, width = 16, height = 8)
g

Figure 4.23: Interpolated satellite tracks at 1 position/3hours of right whales tagged in the Auckland Islands in 2020 and 2021.

g <- gbase + 
  geom_sf(data = nz_fit_ssm_df, aes(col = mode, geometry = geometry), size = 0.8) +
  geom_sf(data = nz_fit_ssm_df[nz_fit_ssm_df$mode == "ARS",], aes(col = mode, geometry = geometry), size = 0.8) +
  theme(legend.position = "right") +
  facet_wrap(~time_step)

ggsave(g, file = "./Outputs/nz_fit_ssm_tracks_ars.png", dpi = 300, width = 16, height = 8)
g

Figure 4.24: Movement mode of right whales tagged in the Auckland Islands in 2020 and 2021.